library(reshape2)
library(ggplot2)
library(plotly)
library(plyr)
library(Rfast)
library(data.table)
library(knitr)
#library(rstudioapi)
#library(here)
library(viridis)# Get path current script is in
#source_path = file.path(here(), 'simulation', 'final_sim', fsep=.Platform$file.sep)
#source_path = getSourceEditorContext()$path
# For knitting
source_path = '/Volumes/MPRG-Neurocode/Users/christoph/pedlr/code/simulation/'
# Source functions required for this script
source(file.path(source_path, 'Beta_pseudo_sim.R', fsep = .Platform$file.sep))
source(file.path(source_path, 'Gaussian_pseudo_sim.R', fsep = .Platform$file.sep))
source(file.path(source_path, 'Create_design.R', fsep = .Platform$file.sep))
source(file.path(source_path, 'Create_miniblock.R', fsep = .Platform$file.sep))
source(file.path(source_path, 'PE_dep_LR_choice_model.R', fsep = .Platform$file.sep))
source(file.path(source_path, 'PE_dep_LR_choice_model_alt.R', fsep = .Platform$file.sep))
source(file.path(source_path, 'PE_dep_LR_choice_model_interdep.R', fsep = .Platform$file.sep))
source(file.path(source_path, 'Softmax_choice.R', fsep = .Platform$file.sep))
source(file.path(source_path, 'PE_dep_LR_single_model.R', fsep = .Platform$file.sep))# Reward lower and upper boundary
reward_space_lb = 1
reward_space_ub = 100
# Set up space for plotting
reward_space_length = 1000
reward_space = seq(reward_space_lb,
reward_space_ub,
length=reward_space_length)# Parameters
gaussian_mean = (100 * 1/6) * 3
gaussian_sd = (100 * 1/6) / 3
# Distribution
gaussian_densities = dnorm(reward_space, gaussian_mean, gaussian_sd)
# Normalize density
gaussian_densities = scale(gaussian_densities) - min(scale(gaussian_densities))# Set up "skewdness" parameter for easier change of slope
beta_skewedness = 5
# Parameters for lower end beta and upper end beta
beta_a_le = beta_skewedness/3
beta_b_le = 2*beta_skewedness/3
beta_a_ue = 2*beta_skewedness/3
beta_b_ue = beta_skewedness/3
# Create distributions
beta_densities_le = dbeta(seq(0,1,length=reward_space_length), beta_a_le,beta_b_le)
beta_densities_ue = dbeta(seq(0,1,length=reward_space_length), beta_a_ue,beta_b_ue)
# Normalize densities
beta_densities_le = scale(beta_densities_le) - min(scale(beta_densities_le))
beta_densities_ue = scale(beta_densities_ue) - min(scale(beta_densities_le))# Form data frame from density values
data_densities = data.frame(reward_space,
gaussian_densities,
beta_densities_le,
beta_densities_ue)
# Bring data frame into long format
data_densities = melt(data_densities, id.vars='reward_space')
# Disribution plot
p_dist = ggplot(data=data_densities, aes(x=reward_space, y=value, group=as.factor(variable),
color=as.factor(variable))) +
geom_line(size=1) +
scale_color_manual(values = c(color$gaussian, color$gain, color$loss)) +
scale_x_continuous(limits = c(1,100)) +
#stat_summary(fun.data = 'mean', geom = 'vline') +
geom_vline(xintercept=c(1/3, 1/2, 2/3)*100, linetype='dashed', alpha=0.5, size=0.1) +
labs(title='Density functions of used distributions',
x='Outcome',
y='Normalized density',
color='Distributions') +
scale_y_continuous(limits=c(0,4.5)) +
theme(plot.title=element_text(size=15, face='bold', hjust=0.5),
axis.title=element_text(size=12))
p_distfile = '/Volumes/MPRG-Neurocode/Users/christoph/PE_dependent_LR/simulation/figure/distributions.pdf'
ggsave(file, plot = p_dist, height=11, width=20, unit='cm')# Disribution plot
p_dist = ggplot(data=data_densities, aes(x=reward_space, y=value, group=as.factor(variable),
color=as.factor(variable),
alpha=as.factor(variable))) +
geom_line(size=1) +
scale_color_brewer(palette='Accent') +
scale_x_continuous(limits = c(1,100)) +
#stat_summary(fun.data = 'mean', geom = 'vline') +
geom_vline(xintercept=c(1/3, 1/2, 2/3)*100, linetype='dashed', alpha=0.5, size=0.1) +
labs(title='Density functions of used distributions',
x='Outcome',
y='Normalized density',
color='Distributions') +
scale_y_continuous(limits=c(0,4.5)) +
theme(plot.title=element_text(size=15, face='bold', hjust=0.5),
axis.title=element_text(size=12))
# Different plots for LIFE academy talk
p_dist_a = p_dist + scale_alpha_manual(values=c(0.3,1,0.3))
file = file.path(source_path, 'figures', 'distributions_a.pdf', fsep = .Platform$file.sep)
ggsave(file, plot = p_dist_a, height=11, width=20, unit='cm')
p_dist_b = p_dist + scale_alpha_manual(values=c(1,0.3,0.3))
file = file.path(source_path, 'figures', 'distributions_b.pdf', fsep = .Platform$file.sep)
ggsave(file, plot = p_dist_b, height=11, width=20, unit='cm')
p_dist_c = p_dist + scale_alpha_manual(values=c(0.3,0.3,1))
file = file.path(source_path, 'figures', 'distributions_c.pdf', fsep = .Platform$file.sep)
ggsave(file, plot = p_dist_c, height=11, width=20, unit='cm')
p_dist_ab = p_dist + scale_alpha_manual(values=c(1,1,0.3))
file = file.path(source_path, 'figures', 'distributions_ab.pdf', fsep = .Platform$file.sep)
ggsave(file, plot = p_dist_ab, height=11, width=20, unit='cm')
p_dist_bc = p_dist + scale_alpha_manual(values=c(1,0.3,1))
file = file.path(source_path, 'figures', 'distributions_bc.pdf', fsep = .Platform$file.sep)
ggsave(file, plot = p_dist_bc, height=11, width=20, unit='cm')# Define parameters
n_samples = 500
data_sampling = data.frame(c(1:n_samples))
colnames(data_sampling) = 'Sample'
# Fill data frame with samples
data_sampling$stim_1 = Beta_pseudo_sim(n_samples, beta_a_le, beta_b_le, 'b_le',
reward_space_lb, reward_space_ub)$outcome
data_sampling$stim_2 = Gaussian_pseudo_sim(n_samples, gaussian_mean, gaussian_sd, 'g_an',
reward_space_lb, reward_space_ub)$outcome
data_sampling$stim_3 = Beta_pseudo_sim(n_samples, beta_a_ue, beta_b_ue, 'b_ue',
reward_space_lb, reward_space_ub)$outcome
data_sampling = melt(data_sampling, id.vars='Sample')
# Plot histograms
p_sampling = ggplot(data=data_sampling, aes(x=value, fill=variable)) +
scale_fill_manual(values=c(color$gain, color$gaussian, color$loss)) +
geom_histogram(binwidth = 1) +
facet_wrap(~variable)
ggplotly(p_sampling)# Cut down data frame for plotting
data_plot = data_model_2b[,1:11]
data_plot = melt(data_plot, id.vars=c('sub_id',
'trial',
'stim_1',
'stim_2',
'reward_stim_1',
'reward_stim_2',
'choice',
'comp_number'))
# Average over all subjects
data_plot = ddply(data_plot,
.(trial, variable),
summarize,
mean_value = mean(value),
sd_value = sd (value))
# Plot model values
p_choice = ggplot(data=data_plot, aes(x=trial, y=mean_value, group=variable, color=variable)) +
geom_hline(yintercept=c(1/3*100, 50, 2/3*100), linetype='dashed') +
geom_hline(yintercept=c(mean(subset(data_plot, variable == 'v_stim_1')$mean_value),
mean(subset(data_plot, variable == 'v_stim_2')$mean_value),
mean(subset(data_plot, variable == 'v_stim_3')$mean_value)),
color = c(color$gain, color$gaussian, color$loss)) +
geom_line(size=1, alpha=0.8) +
scale_color_manual(values = c(color$gain, color$gaussian, color$loss)) +
labs(title=paste('Choice model performance (α0=', alpha0, ', α1=', alpha1, ', τ=', temperature, ')', sep=''),
x='Trial',
y='Estimated outcome',
color='Stimulus')
ggplotly(p_choice)data_plot = data_model_2b[,1:11]
data_plot = melt(data_plot, id.vars=c('sub_id',
'trial',
'stim_1',
'stim_2',
'reward_stim_1',
'reward_stim_2',
'choice',
'comp_number'))
# Select only choices of stim 1 to calculate mean value (leaving out plateaus)
data_1 = subset(data_plot, choice == 1 & variable == 'v_stim_1')
# Calculate mean value of dist 1 WHEN CHOSEN
mean_1 = mean(data_1$value, na.rm = TRUE)
# Same for 2
data_2 = subset(data_plot, choice == 2 & variable == 'v_stim_2')
mean_2 = mean(data_2$value, na.rm = TRUE)
# Same for 3
data_3 = subset(data_plot, choice == 3 & variable == 'v_stim_3')
mean_3 = mean(data_3$value, na.rm = TRUE)
# Average over all subjects
data_plot = ddply(data_plot,
.(trial, variable),
summarize,
mean_value = mean(value),
sd_value = sd (value))
# Plot model values
p = ggplot(data=data_plot, aes(x=trial, y=mean_value, group=variable, color=variable)) +
geom_hline(yintercept=c(1/3*100, 50, 2/3*100), linetype='dashed') +
geom_hline(yintercept=c(mean_1,
mean_2,
mean_3),
color=c(color$gain, color$gaussian, color$loss)) +
geom_line(size=1, alpha=0.8) +
scale_color_manual(values = c(color$gain, color$gaussian, color$loss)) +
labs(title=paste('Choice model performance (α0=', alpha0, ', α1=', alpha1, ', τ=', temperature, ')', sep=''),
x='Trial',
y='Estimated outcome',
color='Stimulus')
p = ggplotly(p)
p# Cut down data frame for plotting
data_plot = data_model_2b[,1:11]
data_plot = melt(data_plot, id.vars=c('sub_id',
'trial',
'stim_1',
'stim_2',
'reward_stim_1',
'reward_stim_2',
'choice',
'comp_number'))
# Select single subject
#data_plot = subset(data_plot, sub_id == 1)
# Plot model values
p_choice = ggplot(data=data_plot, aes(x=trial, y=value, group=variable, color=variable)) +
theme_bw() +
geom_hline(yintercept=c(1/3 *100, 50, 2/3*100), linetype='dashed',
size=0.1) +
# geom_hline(yintercept=c(mean(subset(data_plot, variable == 'v_stim_1')$mean_value),
# mean(subset(data_plot, variable == 'v_stim_2')$mean_value),
# mean(subset(data_plot, variable == 'v_stim_3')$mean_value)),
# color = c(color$gain, color$gaussian, color$loss)) +
geom_line(size=0.1, alpha=0.8) +
scale_color_manual(values=c(color$gain, color$gaussian, color$loss)) +
labs(title=paste('Choice model performance (alpha0=', alpha0, ', alpha1=', alpha1, ', temp=', temperature, ')', sep=''),
x='Trial',
y='Estimated outcome',
color='Stimulus') +
facet_wrap( ~ sub_id, ncol = 8) +
theme(legend.position = 'none',
axis.text = element_text(size=5),
strip.text = element_text(size=5),
panel.grid = element_blank())
ggplotly(p_choice)file = '/Volumes/MPRG-Neurocode/Users/christoph/pedlr/code/simulation/figures/choice_bc.pdf'
ggsave(file, plot = p_choice, height=11, width=20, unit='cm')# Set up matrix to store data in (to append subjects)
data_model_single_gain = data.frame(matrix(0,0,6))
colnames(data_model_single_gain)= c('sub_id', 'trial', 'reward', 'value', 'PE', 'fPE')
# Loop over each subject
for(sub_count in subjects){
# Simulate data for subject
design = Create_design(25,
beta_a_le, beta_b_le,
gaussian_mean, gaussian_sd,
beta_a_ue, beta_b_ue,
two_betas = TRUE,
reward_space_lb, reward_space_ub)
# Extract gain beta rewards
reward_index_left = which(design$option_left == 1)
# Eliminate forced choice trials from left stimulus
reward_index_right = which(design$option_right == 1)
reward_index_right = reward_index_right[-which(reward_index_right %in% reward_index_left)]
reward = c(design$reward_stim_1[reward_index_left], design$reward_stim_2[reward_index_right])
# Fill design with rewards stemming only from gain beta sampling
design = data.frame(matrix(NA, length(reward), 2))
colnames(design) = c('trial', 'reward')
design$trial = c(1:nrow(design))
design$reward = reward
# Create subject specific 'results' data frame to append subjects
results = data.frame(matrix(NA, nrow(design), 6))
colnames(results) = c('sub_id', 'trial', 'reward', 'value', 'PE', 'fPE')
results[,1] = sub_count
results[,2:3] = design
# Run model over simulated data
sim = PE_dep_LR_single_model(design, alpha0, alpha1, temperature)
# Save model values, PE and fPE for ech trial and subject into matrix
results[,4] = sim$values$value[-length(sim$values$value)]
results[,5] = sim$PE$pe
results[,6] = sim$fPE$fpe
# Append all subjects
data_model_single_gain = rbind(data_model_single_gain, results)
}
# Sort after participants
data_model_single_gain = data_model_single_gain[order(data_model_single_gain$sub_id),]# Prepare data frame for plotting
data_plot = data_model_single_gain[1:4]
data_plot = melt(data_plot, id.vars=c('sub_id',
'trial',
'reward'))
# Average over all subjects
data_plot = ddply(data_plot,
.(trial, variable),
summarize,
mean_value = mean(value),
sd_value = sd (value))
# Plot model values
p_choice = ggplot(data=data_plot, aes(x=trial, y=mean_value, group=variable)) +
geom_line(size=1, alpha=0.8,
color=color$gain) +
geom_hline(yintercept=1/3*100, linetype='dashed') +
geom_hline(aes(yintercept=mean(mean_value)),
color=color$gain) +
labs(title=paste('Single model performance (α0=', alpha0, ', α1=', alpha1, ', τ=', temperature, ')', sep=''),
x='Trial',
y='Estimated outcome',
color='Stimulus') +
theme(legend.position = 'none')
ggplotly(p_choice)# Set up matrix to store data in (to append subjects)
data_model_single_loss = data.frame(matrix(0,0,6))
colnames(data_model_single_loss)= c('sub_id', 'trial', 'reward', 'value', 'PE', 'fPE')
# Loop over each subject
for(sub_count in subjects){
# Simulate data for subject
design = Create_design(25,
beta_a_le, beta_b_le,
gaussian_mean, gaussian_sd,
beta_a_ue, beta_b_ue,
two_betas = TRUE,
reward_space_lb, reward_space_ub)
# Extract gain beta rewards
reward_index_left = which(design$option_left == 3)
# Eliminate forced choice trials from left stimulus
reward_index_right = which(design$option_right == 3)
reward_index_right = reward_index_right[-which(reward_index_right %in% reward_index_left)]
reward = c(design$reward_stim_1[reward_index_left], design$reward_stim_2[reward_index_right])
# Fill design with rewards stemming only from gain beta sampling
design = data.frame(matrix(NA, length(reward), 2))
colnames(design) = c('trial', 'reward')
design$trial = c(1:nrow(design))
design$reward = reward
# Create subject specific 'results' data frame to append subjects
results = data.frame(matrix(NA, nrow(design), 6))
colnames(results) = c('sub_id', 'trial', 'reward', 'value', 'PE', 'fPE')
results[,1] = sub_count
results[,2:3] = design
# Run model over simulated data
sim = PE_dep_LR_single_model(design, alpha0, alpha1, temperature)
# Save model values, PE and fPE for ech trial and subject into matrix
results[,4] = sim$values$value[-length(sim$values$value)]
results[,5] = sim$PE$pe
results[,6] = sim$fPE$fpe
# Append all subjects
data_model_single_loss = rbind(data_model_single_loss, results)
}
# Sort after participants
data_model_single_loss = data_model_single_loss[order(data_model_single_loss$sub_id),]# Prepare data frame for plotting
data_plot = data_model_single_loss[1:4]
data_plot = melt(data_plot, id.vars=c('sub_id',
'trial',
'reward'))
# Average over all subjects
data_plot = ddply(data_plot,
.(trial, variable),
summarize,
mean_value = mean(value),
sd_value = sd (value))
# Plot model values
p_choice = ggplot(data=data_plot, aes(x=trial, y=mean_value, group=variable, color=variable)) +
geom_line(size=1, alpha=0.8,
color=color$loss) +
geom_hline(yintercept=2/3*100, linetype='dashed') +
geom_hline(aes(yintercept=mean(mean_value)),
color=color$loss) +
labs(title=paste('Single model performance (α0=', alpha0, ', α1=', alpha1, ', τ=', temperature, ')', sep=''),
x='Trial',
y='Estimated outcome',
color='Stimulus') +
theme(legend.position = 'none')
ggplotly(p_choice)Alternative model which can switch between greedy and softmax (and can use a variable reward space).
# Set up matrix to store data in (to append subjects)
data_model = data.frame(matrix(0,0,17))
# Loop over each subject
for(sub_count in subjects){
set.seed(sub_count)
# Simulate data for subject
design = Create_design(25,
beta_a_le, beta_b_le,
gaussian_mean, gaussian_sd,
beta_a_ue, beta_b_ue,
two_betas = TRUE,
reward_space_lb, reward_space_ub)
# Create matrix to store data for each subject
results = data.frame(matrix(0,nrow(design),ncol(data_model)))
results[,1] = subjects[sub_count]
results[,2] = c(1:nrow(design))
results[,3] = design$option_left
results[,4] = design$option_right
results[,5] = design$reward_stim_1
results[,6] = design$reward_stim_2
results[,7] = design$comp_number
# Run model over simulated data
sim = PE_dep_LR_choice_model_alt(design, alpha0, alpha1, temperature, reward_space_ub, 'softmax')
# Save model values, PE and fPE for ech trial and subject into matrix
results[,8] = sim$choices
results[,9:11] = sim$values
results[,12:14] = sim$PE
results[,15:17] = sim$fPE
# Append all subjects
data_model = rbind(data_model, results)
}
# Convert to data frame (How our data will look like in the study)
data_model = data.frame(data_model)
colnames(data_model) = c('sub_id', 'trial', 'stim_1', 'stim_2',
'reward_stim_1', 'reward_stim_2', 'comp_number', 'choice',
'v_stim_1', 'v_stim_2', 'v_stim_3',
'pe_stim_1', 'pe_stim_2', 'pe_stim_3',
'fpe_stim_1', 'fpe_stim_2', 'fpe_stim_3')
# Sort after participants
data_model = data_model[order(data_model$sub_id),]data_plot = data_model[,1:11]
data_plot = melt(data_plot, id.vars=c('sub_id',
'trial',
'stim_1',
'stim_2',
'reward_stim_1',
'reward_stim_2',
'choice',
'comp_number'))
# Select only choices of stim 1 to calculate mean value (leaving out plateaus)
data_1 = subset(data_plot, choice == 1 & variable == 'v_stim_1')
# Calculate mean value of dist 1 WHEN CHOSEN
mean_1 = mean(data_1$value, na.rm = TRUE)
# Same for 2
data_2 = subset(data_plot, choice == 2 & variable == 'v_stim_2')
mean_2 = mean(data_2$value, na.rm = TRUE)
# Same for 3
data_3 = subset(data_plot, choice == 3 & variable == 'v_stim_3')
mean_3 = mean(data_3$value, na.rm = TRUE)
# Average over all subjects
data_plot = ddply(data_plot,
.(trial, variable),
summarize,
mean_value = mean(value),
sd_value = sd (value))
# Plot model values
p = ggplot(data=data_plot, aes(x=trial, y=mean_value, group=variable, color=variable)) +
geom_hline(yintercept=c(1/3*100, 50, 2/3*100), linetype='dashed') +
geom_hline(yintercept=c(mean_1,
mean_2,
mean_3),
color=c(color$gain, color$gaussian, color$loss)) +
geom_line(size=1, alpha=0.8) +
scale_color_manual(values = c(color$gain, color$gaussian, color$loss)) +
labs(title=paste('Choice model performance (α0=', alpha0, ', α1=', alpha1, ', τ=', temperature, ')', sep=''),
x='Trial',
y='Estimated outcome',
color='Stimulus')
p = ggplotly(p)
p# Set up matrix to store data in (to append subjects)
data_model = data.frame(matrix(0,0,17))
# Loop over each subject
for(sub_count in subjects){
set.seed(sub_count)
# Simulate data for subject
design = Create_design(25,
beta_a_le, beta_b_le,
gaussian_mean, gaussian_sd,
beta_a_ue, beta_b_ue,
two_betas = TRUE,
reward_space_lb, reward_space_ub)
# Create matrix to store data for each subject
results = data.frame(matrix(0,nrow(design),ncol(data_model)))
results[,1] = subjects[sub_count]
results[,2] = c(1:nrow(design))
results[,3] = design$option_left
results[,4] = design$option_right
results[,5] = design$reward_stim_1
results[,6] = design$reward_stim_2
results[,7] = design$comp_number
# Run model over simulated data
sim = PE_dep_LR_choice_model_alt(design, alpha0, alpha1, temperature, reward_space_ub, 'greedy')
# Save model values, PE and fPE for ech trial and subject into matrix
results[,8] = sim$choices
results[,9:11] = sim$values
results[,12:14] = sim$PE
results[,15:17] = sim$fPE
# Append all subjects
data_model = rbind(data_model, results)
}
# Convert to data frame (How our data will look like in the study)
data_model = data.frame(data_model)
colnames(data_model) = c('sub_id', 'trial', 'stim_1', 'stim_2',
'reward_stim_1', 'reward_stim_2', 'comp_number', 'choice',
'v_stim_1', 'v_stim_2', 'v_stim_3',
'pe_stim_1', 'pe_stim_2', 'pe_stim_3',
'fpe_stim_1', 'fpe_stim_2', 'fpe_stim_3')
# Sort after participants
data_model = data_model[order(data_model$sub_id),]data_plot = data_model[,1:11]
data_plot = melt(data_plot, id.vars=c('sub_id',
'trial',
'stim_1',
'stim_2',
'reward_stim_1',
'reward_stim_2',
'choice',
'comp_number'))
# Select only choices of stim 1 to calculate mean value (leaving out plateaus)
data_1 = subset(data_plot, choice == 1 & variable == 'v_stim_1')
# Calculate mean value of dist 1 WHEN CHOSEN
mean_1 = mean(data_1$value, na.rm = TRUE)
# Same for 2
data_2 = subset(data_plot, choice == 2 & variable == 'v_stim_2')
mean_2 = mean(data_2$value, na.rm = TRUE)
# Same for 3
data_3 = subset(data_plot, choice == 3 & variable == 'v_stim_3')
mean_3 = mean(data_3$value, na.rm = TRUE)
# Average over all subjects
data_plot = ddply(data_plot,
.(trial, variable),
summarize,
mean_value = mean(value),
sd_value = sd (value))
# Plot model values
p = ggplot(data=data_plot, aes(x=trial, y=mean_value, group=variable, color=variable)) +
geom_hline(yintercept=c(1/3*100, 50, 2/3*100), linetype='dashed') +
geom_hline(yintercept=c(mean_1,
mean_2,
mean_3),
color=c(color$gain, color$gaussian, color$loss)) +
geom_line(size=1, alpha=0.8) +
scale_color_manual(values = c(color$gain, color$gaussian, color$loss)) +
labs(title=paste('Choice model performance (α0=', alpha0, ', α1=', alpha1, ', τ=', temperature, ')', sep=''),
x='Trial',
y='Estimated outcome',
color='Stimulus')
p = ggplotly(p)
pAdditional parameter handely contribution of different LR components. For now we fix the parameter to 0.5 (equal influence).
# Set up matrix to store data in (to append subjects)
data_model = data.frame(matrix(0,0,17))
# Loop over each subject
for(sub_count in subjects){
set.seed(sub_count)
# Simulate data for subject
design = Create_design(25,
beta_a_le, beta_b_le,
gaussian_mean, gaussian_sd,
beta_a_ue, beta_b_ue,
two_betas = TRUE,
reward_space_lb, reward_space_ub)
# Create matrix to store data for each subject
results = data.frame(matrix(0,nrow(design),ncol(data_model)))
results[,1] = subjects[sub_count]
results[,2] = c(1:nrow(design))
results[,3] = design$option_left
results[,4] = design$option_right
results[,5] = design$reward_stim_1
results[,6] = design$reward_stim_2
results[,7] = design$comp_number
# Run model over simulated data
sim = PE_dep_LR_choice_model_interdep(design,
alpha0,
alpha1,
interdep,
temperature,
reward_space_ub,
'softmax')
# Save model values, PE and fPE for ech trial and subject into matrix
results[,8] = sim$choices
results[,9:11] = sim$values
results[,12:14] = sim$PE
results[,15:17] = sim$fPE
# Append all subjects
data_model = rbind(data_model, results)
}
# Convert to data frame (How our data will look like in the study)
data_model = data.frame(data_model)
colnames(data_model) = c('sub_id', 'trial', 'stim_1', 'stim_2',
'reward_stim_1', 'reward_stim_2', 'comp_number', 'choice',
'v_stim_1', 'v_stim_2', 'v_stim_3',
'pe_stim_1', 'pe_stim_2', 'pe_stim_3',
'fpe_stim_1', 'fpe_stim_2', 'fpe_stim_3')
# Sort after participants
data_model = data_model[order(data_model$sub_id),]data_plot = data_model[,1:11]
data_plot = melt(data_plot, id.vars=c('sub_id',
'trial',
'stim_1',
'stim_2',
'reward_stim_1',
'reward_stim_2',
'choice',
'comp_number'))
# Select only choices of stim 1 to calculate mean value (leaving out plateaus)
data_1 = subset(data_plot, choice == 1 & variable == 'v_stim_1')
# Calculate mean value of dist 1 WHEN CHOSEN
mean_1 = mean(data_1$value, na.rm = TRUE)
# Same for 2
data_2 = subset(data_plot, choice == 2 & variable == 'v_stim_2')
mean_2 = mean(data_2$value, na.rm = TRUE)
# Same for 3
data_3 = subset(data_plot, choice == 3 & variable == 'v_stim_3')
mean_3 = mean(data_3$value, na.rm = TRUE)
# Average over all subjects
data_plot = ddply(data_plot,
.(trial, variable),
summarize,
mean_value = mean(value),
sd_value = sd (value))
# Plot model values
p = ggplot(data=data_plot, aes(x=trial, y=mean_value, group=variable, color=variable)) +
geom_hline(yintercept=c(1/3*100, 50, 2/3*100), linetype='dashed') +
geom_hline(yintercept=c(mean_1,
mean_2,
mean_3),
color=c(color$gain, color$gaussian, color$loss)) +
geom_line(size=1, alpha=0.8) +
scale_color_manual(values = c(color$gain, color$gaussian, color$loss)) +
labs(title=paste('Choice model performance (α0=', alpha0, ', α1=', alpha1, ', τ=', temperature, ')', sep=''),
x='Trial',
y='Estimated outcome',
color='Stimulus')
p = ggplotly(p)
p# Set up matrix to store data in (to append subjects)
data_model = matrix(0,0,17)
# Loop over each subject
for(sub_count in subjects){
set.seed(sub_count)
# Simulate data for subject
design = Create_design(25,
beta_a_le, beta_b_le,
gaussian_mean, gaussian_sd,
beta_a_ue, beta_b_ue,
two_betas = TRUE,
reward_space_lb, reward_space_ub)
# Switch positions of distributions (gain beta on UE, loss beta on LE)
# Gain beta
reward_index_left = which(design$option_left == 1)
design$reward_stim_1[reward_index_left] = design$reward_stim_1[reward_index_left] + 33
design$reward_stim_1[which(design$reward_stim_1 >= 100)] = 100
reward_index_right = which(design$option_right == 1)
design$reward_stim_2[reward_index_right] = design$reward_stim_2[reward_index_right] + 33
design$reward_stim_2[which(design$reward_stim_2 >= 100)] = 100
# Loss beta
reward_index_left = which(design$option_left == 3)
design$reward_stim_1[reward_index_left] = design$reward_stim_1[reward_index_left] - 33
design$reward_stim_1[which(design$reward_stim_1 <= 1)] = 1
reward_index_right = which(design$option_right == 3)
design$reward_stim_2[reward_index_right] = design$reward_stim_2[reward_index_right] - 33
design$reward_stim_2[which(design$reward_stim_2 <= 1)] = 1
# Create matrix to store data for each subject
results = matrix(0,nrow(design),ncol(data_model))
results[,1] = subjects[sub_count]
results[,2] = c(1:nrow(design))
results[,3] = design$option_left
results[,4] = design$option_right
results[,5] = design$reward_stim_1
results[,6] = design$reward_stim_2
results[,7] = design$comp_number
# Run model over simulated data
sim = PE_dep_LR_choice_model(design, alpha0, alpha1, temperature)
# Save model values, PE and fPE for ech trial and subject into matrix
results[,8] = sim$choice$choice
results[,9:11] = as.matrix(sim$values)
results[,12:14] = as.matrix(sim$PE)
results[,15:17] = as.matrix(sim$fPE)
# Append all subjects
data_model = rbind(data_model, results)
}
# Convert to data frame (How our data will look like in the study)
data_model = data.frame(data_model)
colnames(data_model) = c('sub_id', 'trial', 'stim_1', 'stim_2',
'reward_stim_1', 'reward_stim_2', 'comp_number', 'choice',
'v_stim_1', 'v_stim_2', 'v_stim_3',
'pe_stim_1', 'pe_stim_2', 'pe_stim_3',
'fpe_stim_1', 'fpe_stim_2', 'fpe_stim_3')
# Sort after participants
data_model = data_model[order(data_model$sub_id),]gain_samples = c(subset(data_model, 'stim_1' != 'stim_2')$reward_stim_1[subset(data_model, 'stim_1' != 'stim_2')$stim_1 == 1],
subset(data_model, 'stim_1' != 'stim_2')$reward_stim_2[subset(data_model, 'stim_1' != 'stim_2')$stim_2 == 1],
subset(data_model, 'stim_1' == 'stim_2')$reward_stim_2[subset(data_model, 'stim_1' == 'stim_2')$stim_1 == 1])
loss_samples = c(subset(data_model, 'stim_1' != 'stim_2')$reward_stim_1[subset(data_model, 'stim_1' != 'stim_2')$stim_1 == 3],
subset(data_model, 'stim_1' != 'stim_2')$reward_stim_2[subset(data_model, 'stim_1' != 'stim_2')$stim_2 == 3],
subset(data_model, 'stim_1' == 'stim_2')$reward_stim_2[subset(data_model, 'stim_1' == 'stim_2')$stim_1 == 3])
gauss_samples = c(subset(data_model, 'stim_1' != 'stim_2')$reward_stim_1[subset(data_model, 'stim_1' != 'stim_2')$stim_1 == 2],
subset(data_model, 'stim_1' != 'stim_2')$reward_stim_2[subset(data_model, 'stim_1' != 'stim_2')$stim_2 == 2],
subset(data_model, 'stim_1' == 'stim_2')$reward_stim_2[subset(data_model, 'stim_1' == 'stim_2')$stim_1 == 2])
data_sampling = data.frame(matrix(NA, length(gain_samples), 3))
colnames(data_sampling) = c('stim_1', 'stim_2', 'stim_3')
data_sampling$stim_1 = gain_samples
data_sampling$stim_2 = gauss_samples
data_sampling$stim_3 = loss_samples
data_sampling = melt(data_sampling)
# Plot histograms
p_sampling = ggplot(data=data_sampling, aes(x=value, fill=variable)) +
scale_fill_manual(values = c(color$gain, color$gaussian, color$loss)) +
geom_histogram(binwidth = 1) +
facet_wrap(~variable)
ggplotly(p_sampling)The gain and loss beta distributions are now on the swapped ends of the reward space. If this would switch the bias we would know its due to the position and maybe the undersampling.
data_plot = data_model[,1:11]
data_plot = melt(data_plot, id.vars=c('sub_id',
'trial',
'stim_1',
'stim_2',
'reward_stim_1',
'reward_stim_2',
'choice',
'comp_number'))
# Select only choices of stim 1 to calculate mean value (leaving out plateaus)
data_1 = subset(data_plot, choice == 1 & variable == 'v_stim_1')
# Calculate mean value of dist 1 WHEN CHOSEN
mean_1 = mean(data_1$value, na.rm = TRUE)
# Same for 2
data_2 = subset(data_plot, choice == 2 & variable == 'v_stim_2')
mean_2 = mean(data_2$value, na.rm = TRUE)
# Same for 3
data_3 = subset(data_plot, choice == 3 & variable == 'v_stim_3')
mean_3 = mean(data_3$value, na.rm = TRUE)
# Average over all subjects
data_plot = ddply(data_plot,
.(trial, variable),
summarize,
mean_value = mean(value),
sd_value = sd (value))
# Plot model values
p = ggplot(data=data_plot, aes(x=trial, y=mean_value, group=variable, color=variable)) +
geom_hline(yintercept=c(1/3*100, 50, 2/3*100), linetype='dashed') +
geom_hline(yintercept=c(mean_1,
mean_2,
mean_3),
color=c(color$gain, color$gaussian, color$loss)) +
geom_line(size=1, alpha=0.8) +
scale_color_manual(values = c(color$gain, color$gaussian, color$loss)) +
labs(title=paste('Choice model performance (α0=', alpha0, ', α1=', alpha1, ', τ=', temperature, ')', sep=''),
x='Trial',
y='Estimated outcome',
color='Stimulus')
p = ggplotly(p)
pThe stimuli changed position but the bias pattern stayed the same. The bias does not swap for the loss beta, even if it is posted to the lower end. Same for the gain beta.
# Set up matrix to store data in (to append subjects)
data_model = matrix(0,0,17)
# Loop over each subject
for(sub_count in subjects){
set.seed(sub_count)
# Simulate data for subject
design = Create_design(25,
beta_a_le, beta_b_le,
gaussian_mean, gaussian_sd,
beta_a_ue, beta_b_ue,
two_betas = TRUE,
reward_space_lb, reward_space_ub)
# Resample beta distirbutions as Gaussians
# Gain beta
# Find gain beta comparisons
reward_index_left = which(design$option_left == 1)
reward_index_right = which(design$option_right == 1)
# Delete forced choices from comparisons to handle them separately
reward_index_forced = intersect(reward_index_left, reward_index_right)
reward_index_left = reward_index_left[!reward_index_left %in% reward_index_forced]
reward_index_right = reward_index_right[!reward_index_right %in% reward_index_forced]
# Sample as many rewards as needed from LE gaussian
replace = Gaussian_pseudo_sim(length(c(reward_index_left, reward_index_right, reward_index_forced)),
1/3 * 100,
gaussian_sd,
'gauss_ue',
0,
100)
# Replace rewards of gain beta
design$reward_stim_1[reward_index_left] = replace$outcome[1:length(reward_index_left)]
design$reward_stim_2[reward_index_right] = replace$outcome[(length(reward_index_left)+1):(length(reward_index_right)+length(reward_index_left))]
# Forced choices are copied for both stimuli
design$reward_stim_1[reward_index_forced] = replace$outcome[(length(reward_index_right)+
length(reward_index_left)+1):(
length(reward_index_right)+
length(reward_index_left)+
length(reward_index_forced))]
design$reward_stim_2[reward_index_forced] = design$reward_stim_1[reward_index_forced]
# Loss beta
reward_index_left = which(design$option_left == 3)
reward_index_right = which(design$option_right == 3)
# Delete forced choices from comparisons to handle them separately
reward_index_forced = intersect(reward_index_left, reward_index_right)
reward_index_left = reward_index_left[!reward_index_left %in% reward_index_forced]
reward_index_right = reward_index_right[!reward_index_right %in% reward_index_forced]
# Sample as many rewards as needed from LE gaussian
replace = Gaussian_pseudo_sim(length(c(reward_index_left, reward_index_right, reward_index_forced)),
2/3 * 100,
gaussian_sd,
'gauss_ue',
0,
100)
# Replace rewards of gain beta
# Replace rewards of gain beta
design$reward_stim_1[reward_index_left] = replace$outcome[1:length(reward_index_left)]
design$reward_stim_2[reward_index_right] = replace$outcome[(length(reward_index_left)+1):(length(reward_index_right)+length(reward_index_left))]
# Forced choices are copied for both stimuli
design$reward_stim_1[reward_index_forced] = replace$outcome[(length(reward_index_right)+
length(reward_index_left)+1):(
length(reward_index_right)+
length(reward_index_left)+
length(reward_index_forced))]
design$reward_stim_2[reward_index_forced] = design$reward_stim_1[reward_index_forced]
# Create matrix to store data for each subject
results = matrix(0,nrow(design),ncol(data_model))
results[,1] = subjects[sub_count]
results[,2] = c(1:nrow(design))
results[,3] = design$option_left
results[,4] = design$option_right
results[,5] = design$reward_stim_1
results[,6] = design$reward_stim_2
results[,7] = design$comp_number
# Run model over simulated data
sim = PE_dep_LR_choice_model(design, alpha0, alpha1, temperature)
# Save model values, PE and fPE for ech trial and subject into matrix
results[,8] = sim$choice$choice
results[,9:11] = as.matrix(sim$values)
results[,12:14] = as.matrix(sim$PE)
results[,15:17] = as.matrix(sim$fPE)
# Append all subjects
data_model = rbind(data_model, results)
}
# Convert to data frame (How our data will look like in the study)
data_model = data.frame(data_model)
colnames(data_model) = c('sub_id', 'trial', 'stim_1', 'stim_2',
'reward_stim_1', 'reward_stim_2', 'comp_number', 'choice',
'v_stim_1', 'v_stim_2', 'v_stim_3',
'pe_stim_1', 'pe_stim_2', 'pe_stim_3',
'fpe_stim_1', 'fpe_stim_2', 'fpe_stim_3')
# Sort after participants
data_model = data_model[order(data_model$sub_id),]gain_samples = c(subset(data_model, 'stim_1' != 'stim_2')$reward_stim_1[subset(data_model, 'stim_1' != 'stim_2')$stim_1 == 1],
subset(data_model, 'stim_1' != 'stim_2')$reward_stim_2[subset(data_model, 'stim_1' != 'stim_2')$stim_2 == 1],
subset(data_model, 'stim_1' == 'stim_2')$reward_stim_2[subset(data_model, 'stim_1' == 'stim_2')$stim_1 == 1])
loss_samples = c(subset(data_model, 'stim_1' != 'stim_2')$reward_stim_1[subset(data_model, 'stim_1' != 'stim_2')$stim_1 == 3],
subset(data_model, 'stim_1' != 'stim_2')$reward_stim_2[subset(data_model, 'stim_1' != 'stim_2')$stim_2 == 3],
subset(data_model, 'stim_1' == 'stim_2')$reward_stim_2[subset(data_model, 'stim_1' == 'stim_2')$stim_1 == 3])
gauss_samples = c(subset(data_model, 'stim_1' != 'stim_2')$reward_stim_1[subset(data_model, 'stim_1' != 'stim_2')$stim_1 == 2],
subset(data_model, 'stim_1' != 'stim_2')$reward_stim_2[subset(data_model, 'stim_1' != 'stim_2')$stim_2 == 2],
subset(data_model, 'stim_1' == 'stim_2')$reward_stim_2[subset(data_model, 'stim_1' == 'stim_2')$stim_1 == 2])
data_sampling = data.frame(matrix(NA, length(gain_samples), 3))
colnames(data_sampling) = c('stim_1', 'stim_2', 'stim_3')
data_sampling$stim_1 = gain_samples
data_sampling$stim_2 = gauss_samples
data_sampling$stim_3 = loss_samples
data_sampling = melt(data_sampling)
# Plot histograms
p_sampling = ggplot(data=data_sampling, aes(x=value, fill=variable)) +
scale_fill_brewer(palette = 'Accent') +
geom_histogram(binwidth = 1) +
facet_wrap(~variable)
ggplotly(p_sampling)data_plot = data_model[,1:11]
data_plot = melt(data_plot, id.vars=c('sub_id',
'trial',
'stim_1',
'stim_2',
'reward_stim_1',
'reward_stim_2',
'choice',
'comp_number'))
# Select only choices of stim 1 to calculate mean value (leaving out plateaus)
data_1 = subset(data_plot, choice == 1 & variable == 'v_stim_1')
# Calculate mean value of dist 1 WHEN CHOSEN
mean_1 = mean(data_1$value, na.rm = TRUE)
# Same for 2
data_2 = subset(data_plot, choice == 2 & variable == 'v_stim_2')
mean_2 = mean(data_2$value, na.rm = TRUE)
# Same for 3
data_3 = subset(data_plot, choice == 3 & variable == 'v_stim_3')
mean_3 = mean(data_3$value, na.rm = TRUE)
# Average over all subjects
data_plot = ddply(data_plot,
.(trial, variable),
summarize,
mean_value = mean(value),
sd_value = sd (value))
# Plot model values
p = ggplot(data=data_plot, aes(x=trial, y=mean_value, group=variable, color=variable)) +
geom_hline(yintercept=c(1/3*100, 50, 2/3*100), linetype='dashed') +
geom_hline(yintercept=c(mean_1,
mean_2,
mean_3),
color=c(color$gain, color$gaussian, color$loss)) +
geom_line(size=1, alpha=0.8) +
scale_color_manual(values = c(color$gain, color$gaussian, color$loss)) +
labs(title=paste('Choice model performance (α0=', alpha0, ', α1=', alpha1, ', τ=', temperature, ')', sep=''),
x='Trial',
y='Estimated outcome',
color='Stimulus')
p = ggplotly(p)
pBias disappears for Gaussian distributions
Simulate data
Plot results of choice model
data_plot = data_model[,1:11]
data_plot = melt(data_plot, id.vars=c('sub_id',
'trial',
'stim_1',
'stim_2',
'reward_stim_1',
'reward_stim_2',
'choice',
'comp_number'))
# Select only choices of stim 1 to calculate mean value (leaving out plateaus)
data_1 = subset(data_plot, choice == 1 & variable == 'v_stim_1')
# Calculate mean value of dist 1 WHEN CHOSEN
mean_1 = mean(data_1$value, na.rm = TRUE)
# Same for 2
data_2 = subset(data_plot, choice == 2 & variable == 'v_stim_2')
mean_2 = mean(data_2$value, na.rm = TRUE)
# Same for 3
data_3 = subset(data_plot, choice == 3 & variable == 'v_stim_3')
mean_3 = mean(data_3$value, na.rm = TRUE)
# Average over all subjects
data_plot = ddply(data_plot,
.(trial, variable),
summarize,
mean_value = mean(value),
sd_value = sd (value))
# Plot model values
p = ggplot(data=data_plot, aes(x=trial, y=mean_value, group=variable, color=variable)) +
geom_hline(yintercept=c(1/6*100, 2/6*100, 50), linetype='dashed') +
geom_hline(yintercept=c(mean_1,
mean_2,
mean_3),
color=c(color$gaussian, color$gain, color$gaussian)) +
geom_line(size=1, alpha=0.8) +
scale_color_manual(values = c(color$gaussian, color$gain, color$gaussian)) +
labs(title=paste('Choice model performance (α0=', alpha0, ', α1=', alpha1, ', τ=', temperature, ')', sep=''),
x='Trial',
y='Estimated outcome',
color='Stimulus')
p = ggplotly(p)
pSimulate data
Plot results
data_plot = data_model[,1:11]
data_plot = melt(data_plot, id.vars=c('sub_id',
'trial',
'stim_1',
'stim_2',
'reward_stim_1',
'reward_stim_2',
'choice',
'comp_number'))
# Select only choices of stim 1 to calculate mean value (leaving out plateaus)
data_1 = subset(data_plot, choice == 1 & variable == 'v_stim_1')
# Calculate mean value of dist 1 WHEN CHOSEN
mean_1 = mean(data_1$value, na.rm = TRUE)
# Same for 2
data_2 = subset(data_plot, choice == 2 & variable == 'v_stim_2')
mean_2 = mean(data_2$value, na.rm = TRUE)
# Same for 3
data_3 = subset(data_plot, choice == 3 & variable == 'v_stim_3')
mean_3 = mean(data_3$value, na.rm = TRUE)
# Average over all subjects
data_plot = ddply(data_plot,
.(trial, variable),
summarize,
mean_value = mean(value),
sd_value = sd (value))
# Plot model values
p = ggplot(data=data_plot, aes(x=trial, y=mean_value, group=variable, color=variable)) +
geom_hline(yintercept=c(50, 4/6*100, 5/6*100), linetype='dashed') +
geom_hline(yintercept=c(mean_1,
mean_2,
mean_3),
color=c(color$gaussian, color$loss, color$gaussian)) +
geom_line(size=1, alpha=0.8) +
scale_color_manual(values = c(color$gaussian, color$loss, color$gaussian)) +
labs(title=paste('Choice model performance (α0=', alpha0, ', α1=', alpha1, ', τ=', temperature, ')', sep=''),
x='Trial',
y='Estimated outcome',
color='Stimulus')
p = ggplotly(p)
p# Number of samples for mean pe
i = 10
pe_mean_gain = c(1:i)
for(i in 1:10){
# Define simulation parameters
n_subjects = 40
subjects = c(1:n_subjects)
n_miniblocks = 25
# Set up matrix to store data in (to append subjects)
data_model = matrix(0,0,17)
# Loop over each subject
for(sub_count in subjects){
set.seed(sub_count)
# Simulate data for subject
design = Create_design(25,
1/6*100, gaussian_sd,
beta_a_le, beta_b_le,
3/6*100, gaussian_sd,
two_betas = FALSE,
reward_space_lb, reward_space_ub)
# Create matrix to store data for each subject
results = data.frame(matrix(0,nrow(design),ncol(data_model)))
results[,1] = subjects[sub_count]
results[,2] = c(1:nrow(design))
results[,3] = design$option_left
results[,4] = design$option_right
results[,5] = design$reward_stim_1
results[,6] = design$reward_stim_2
results[,7] = design$comp_number
# Run model over simulated data
sim = PE_dep_LR_choice_model_alt(design, alpha0, alpha1, temperature, reward_space_ub, 'softmax')
# Save model values, PE and fPE for ech trial and subject into matrix
results[,8] = sim$choices
results[,9:11] = sim$values
results[,12:14] = sim$PE
results[,15:17] = sim$fPE
# Append all subjects
data_model = rbind(data_model, results)
}
# Convert to data frame (How our data will look like in the study)
data_model = data.frame(data_model)
colnames(data_model) = c('sub_id', 'trial', 'stim_1', 'stim_2',
'reward_stim_1', 'reward_stim_2', 'comp_number', 'choice',
'v_stim_1', 'v_stim_2', 'v_stim_3',
'pe_stim_1', 'pe_stim_2', 'pe_stim_3',
'fpe_stim_1', 'fpe_stim_2', 'fpe_stim_3')
# Sort after participants
data_model = data_model[order(data_model$sub_id),]
# PE analysis
pe_mean_gain[i] = mean(data_model$pe_stim_2[data_model$pe_stim_2 != 0], na.rm = TRUE)
}
pe_mean_gain## [1] -0.5559428 -0.5559428 -0.5559428 -0.5559428 -0.5559428 -0.5559428
## [7] -0.5559428 -0.5559428 -0.5559428 -0.5559428
## [1] -0.5559428
Negative average PE in Gain beta
# Number of samples for mean pe
i = 10
pe_mean_loss = c(1:i)
for(i in 1:10){
# Define simulation parameters
n_subjects = 40
subjects = c(1:n_subjects)
n_miniblocks = 25
# Set up matrix to store data in (to append subjects)
data_model = matrix(0,0,17)
# Loop over each subject
for(sub_count in subjects){
set.seed(sub_count)
# Simulate data for subject
design = Create_design(25,
3/6*100, gaussian_sd,
beta_a_ue, beta_b_ue,
5/6*100, gaussian_sd,
two_betas = FALSE,
reward_space_lb, reward_space_ub)
# Create matrix to store data for each subject
results = data.frame(matrix(0,nrow(design),ncol(data_model)))
results[,1] = subjects[sub_count]
results[,2] = c(1:nrow(design))
results[,3] = design$option_left
results[,4] = design$option_right
results[,5] = design$reward_stim_1
results[,6] = design$reward_stim_2
results[,7] = design$comp_number
# Run model over simulated data
sim = PE_dep_LR_choice_model_alt(design, alpha0, alpha1, temperature, reward_space_ub, 'softmax')
# Save model values, PE and fPE for ech trial and subject into matrix
results[,8] = sim$choices
results[,9:11] = sim$values
results[,12:14] = sim$PE
results[,15:17] = sim$fPE
# Append all subjects
data_model = rbind(data_model, results)
}
# Convert to data frame (How our data will look like in the study)
data_model = data.frame(data_model)
colnames(data_model) = c('sub_id', 'trial', 'stim_1', 'stim_2',
'reward_stim_1', 'reward_stim_2', 'comp_number', 'choice',
'v_stim_1', 'v_stim_2', 'v_stim_3',
'pe_stim_1', 'pe_stim_2', 'pe_stim_3',
'fpe_stim_1', 'fpe_stim_2', 'fpe_stim_3')
# Sort after participants
data_model = data_model[order(data_model$sub_id),]
# PE analysis
pe_mean_loss[i] = mean(data_model$pe_stim_2[data_model$pe_stim_2 != 0], na.rm = TRUE)
}
pe_mean_loss## [1] 0.5552614 0.5552614 0.5552614 0.5552614 0.5552614 0.5552614 0.5552614
## [8] 0.5552614 0.5552614 0.5552614
## [1] 0.5552614
Positive average PE in Loss beta
pe_stim1 = subset(data_model, pe_stim_1 !=0)
update_stim1 = pe_stim1$pe_stim_1 * pe_stim1$fpe_stim_1
mean(update_stim1)## [1] -0.005178617
pe_stim3 = subset(data_model, pe_stim_3 !=0)
update_stim3 = pe_stim3$pe_stim_3 * pe_stim3$fpe_stim_3
mean(update_stim3)## [1] 0.07918103